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Abstract. - We study fixed density sandpiles in which the number of particles transferred to 
a neighbor on relaxing an active site is determined stochastically by a parameter p. Using an 
argument, the critical density at which an active-absorbing transition occurs is found exactly. 
We study the critical behavior numerically and find that the exponents associated with both 
static and time-dependent quantities vary continuously with p. 



An essential goal of the study of nonequilibrium critical phenomena is the identification and 
characterisation of the various universality classes as in the equilibrium case. An important 
prototype of nonequilibrium phase transitions for which extensive research has been done 
in this direction are the ones that occur between active states with sustained activity and 
absorbing states devoid of activity in which the system remains trapped. Typically the critical 
behavior of such systems falls in the directed percolation class but additional symmetries such 
as parity conservation and symmetry among absorbing states are known to change it [1]. 
Recent studies of active-absorbing transition in a class of fixed density sandpile models have 
shown the existence of a new universality class arising due to the conservation of global 
density [2,3]. 

Traditional sandpile models, however, are driven-dissipative systems in which on exceeding 
a local threshold, a site becomes active (or unstable) and relaxes by redistributing its particles 
according to some rules. Most studies have focused on the self-organised critical (SOC) state 
whose behavior depends on the details of the model; for instance, models with deterministic [4] 
or stochastic [5] motion of particles, stochastic thresholds [6] and "sticky grains" due to which 
an active site may or may not relax [7] have different critical exponents. By fixing the total 
density in these systems, the phase transition can be studied and the behavior at the critical 
point is found, under certain conditions [8], to be same as that of the SOC state of the 
corresponding slowly driven sandpile [9,10]. 

A common feature of the sandpiles mentioned above is that the number of particles leaving 
an active site on relaxation is fixed. However, in experimental situations, this need not be the 
case. To this end, we introduce a class of sandpiles in which the number of particles transferred 
to a neighbor is determined stochastically. In the steady state, an active-absorbing transition 
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occurs as the (conserved) total density of the system is tuned. We study the critical behavior 
of an order parameter in the steady state and its relaxation dynamics, and find that the 
associated exponents vary continuously with the stochasticity parameter. 

Model. - Our model is defined on a ring with L sites with a non-negative integer rii of 
particles at each site. We call a site i active if n, is greater than the threshold n c = 1 otherwise 
inactive. A configuration C = {ni,ri2, is updated in continuous time according to the 

following rules. With probability p, one particle hops out of an active site to either left or 
right neighbor with equal probability 

{...,n i _i,n i ,...} ^ {...,m-i + l,m - 1, ...} 

p/2 

{...,rii,n i+ i,...} — ► {...,rii - l,n i+ i + 1, ...} , > 2 , (1) 

whereas with probability q = 1 — p, two particles leave an active site, one to the left and 
another to the right 

{...,Hi-i,rH,ni + i, ...} {...,n,_i + 1,71, - 2,n i+1 + 1, ...} , > 2 . (2) 

Clearly, the total density p = N/L, where N is the total number of particles, is conserved. For 
any < p < 1, the steady state of this system shows a phase transition at the critical density 
p c between an absorbing phase in which each site has either none or one particle and an 
active phase with nonzero density S(p,p) of active sites. Wc find numerically that the critical 
exponents for both the static and the dynamic quantities vary continuously with parameter p. 
Although such nonunivcrsal behavior has been seen in anomalous directed percolation which 
models epidemics with long range infections [11] and in a non-local sandpile [12], the model 
studied here has only on-site interactions. 

Many properties of our system are known exactly in d dimensions for p = and 1. At 
p = 0, two particles leave an active site and move deterministically to each of the neighbors. 
It is known that the SOC state of this model has density p c = (L — l)/L and the perturbations 
about this density relax in typical time L z with dynamic exponent z = 1 [13, 14]. At p = 1, a 
single particle hops out to a neighbor which is determined stochastically. It was shown in [15] 
that in this case, an active-absorbing transition occurs at p c = 1 and due to the absence 
of correlations in the active phase, the activity S(p,l) oc p — p c . Further, a mapping to a 
two-species annihilation problem was used to deduce that at the critical point, the activity 
S(t,p,l) at time t decays as a power law with exponent = 1/4 until a typical time L z 
with z = 2 while in the absorbing phase, it relaxes as a stretched exponential e -1 -*/* ) with 
exponent a = 1/3. In this article, we are interested in the behavior of the model described 
above for < p < 1. 

Critical behavior in the steady state. - We first argue that the critical density p c = 1 
for p > and also give supporting numerical evidence. Since for p > 1, we have at least one 
active site, it follows that p c < 1. At p = 1, consider the rate of change of the probability of 
the configuration Co in which each site has one particle. In the steady state, since the system 
cannot leave this configuration, the probability of configurations entering it must be zero. In 
fact, all the configurations that can reach Co in finite number of time steps must have zero 
probability. For p > 0, any initial configuration can reach Co at large times by filling up 
its empty sites without creating new ones due to the single particle move (eq.Q). Thus, Cq 
is the only allowed configuration at unit density and since the activity S is zero for Co, the 
critical density p c = 1 for Q < p < 1. As shown in fig. ^ for various values of p, at large times, 
the activity S(t,p,p) at p = 1 decays to zero which confirms our assertion that p c = 1. 
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Fig. 1 - Temporal decay of the activity S(t,p,p) at p = 1 with L = 512 for various values of p. 
The activity finally goes to zero after an initial power law regime with p-dependent exponent 9(p), 
9(1) — 1/4. Inset: Data collapse of the activity in ea. HH with z = 5/4 at p = 1/4. 



Having determined the critical density, we are now in a position to study the critical 
behavior of the activity S = ^2 n>2 P( n ) where P(n) is the distribution of the number of 
particles at a site in the steady state. Figure |2] shows our results obtained using Monte Carlo 
simulations for the steady state activity S(p,p) which grows as a power law with deviation 
p — p c from the critical point 

S(p,p)K(p- Pc f , (3) 

where p is an increasing function of p with (3(1) = 1. We have verified that the above scaling 
behavior does not depend on the initial conditions. A useful check on the numerics was 
provided by the following exact relation 

The above equation is just the statement of particle conservation and is valid for all p. At 
p = 1, due to the absence of empty sites in the active phase [15], the normalization condition 
is P(l) + 5 = 1. Using this in eq.(@J, we immediately obtain S oc (p — 1). For p < 1, 
the probability P(0) is nonzero due to eq.J5J) so that in order to determine (3(p), one more 
nontrivial equation is required. The above relation, however, has the interesting consequence 
that both P(0) and 1 — -P(l) also scale with p — p c with exponent j3(p) (see inset of fig. El- 
This scaling behavior can be shown for 1 — P(l) by expanding eq.J3J about the critical density 
and using that (3(p) < 1, and for P(0) the desired result follows on using the normalization. 

To derive eq.(@J), we need to consider the equations obeyed by a connected two-point 
correlation function in the steady state; however, this calculation is cumbersome and lengthy, 
so we do not reproduce the details here and defer it to a future publication [16]. We mention 
that for p < 1, our analysis also shows that the aforementioned correlation function vanishes, 
in the thermodynamic limit, for all but nearest neighbors indicating short range correlations 
in the active phase. Recall that at p = 1 even these correlations are absent [15]. 

To assess the importance of the correlations in the system, we next consider a mean field 
theory which completely ignores the correlations. We start with the exact equations for the 
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Fig. 2 - Activity S(p,p) vs. p - p c for p = l/8(+), l/4(x), 1/2(A), 5/8(o) and 3/4(») with L = 256. 
The curve (p — p c )/p shown with solid line is the exact expression for S at p — 1. Inset: Scaling of 
P(0)(A), and 1 - P(l)(x) with p - p c for p = 1/4 (top) and p = 3/4 (bottom) with L = 2048. 



rate of change of the number distribution P{n, t) at time t. In the steady state, this distribution 
is independent of time and we obtain 



P P{n + l) + qP{n + 2)+p'^P{n-l in ') = P(n)(l - £„,i) + p' ^ P(n,n') t n > 1(5) 

n'>2 n>>2 

qP(2) = p'^P(0,n') , (6) 

n'>2 

where p' = 2 — p and P(n, n') is the joint distribution at two adjacent sites. In the above 
equations, the left hand side represents the elementary moves due to which n particles can be 
obtained at a site, while the right hand side corresponds to the ways by which the particle 
number changes at a site having n particles. 

By approximating the joint distribution function P(n, n') in Eqs.Q and © by the product 
P(n)P(n'), we find that the generating function G(w) = X)n>2 P(n-)w n is given by 

_ ^(i-s-d-.w)) 

(w+ — w)(w — w-) 

where 



l±y/l + 4(l-p)(2-p)S 

W± = 2T2—pjS ' (8) 

with the normalisation G(l) = S. Due to the particle conservation constraint 

P , 1U <K?I (3-2p)S+(2-p)P(l) 

" = P(1) + ^ U=1 = (2-p)(l-g) • (9) 

As expected, this equation coincides with the exact eq.(QJ only at p = 1. In order to solve for S 

as a function of p, we also need to determine P(l). Since G(w) in eq.Q) has a pole at |io_ | < 1 
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Fig. 3 - Scaled relaxation time distribution LP(t, L,p)/p at the critical point with L — 128 for various 
values of p. The lines are for p = 1/4 (top) and p — 3/4 (bottom) and have slope equal to p/(l +p). 
Inset: Scaling collapse of the distribution P(t,L,p) in eq. iffill with z = 3/2 at p = 1/2. 



which corresponds to exponentially growing P{n) with n, we avoid this unphysical solution 
by demanding that tu_ is a root of the numerator as well and obtain P(l) = (1 — S)/(l — tu_). 
Using this expression for P(l) in eq.© and that the activity S is zero at the critical density, 
we find p c = 1/(2 — p) and a Taylor series expansion about this p c yields S « p c {p — p c ). 

Since these mean field predictions do not match the numerical results (except at p = 1), we 
include some correlations by implementing an improved mean field theory which was successful 
in determining the fundamental diagram of a traffic model with parallel update [17]. We first 
note that for p < 1, a site with two particles can be emptied with probability 1 — p but this 
move fills the adjacent (possibly empty) sites so that in the active phase, it is not possible to 
create a pair (and hence clusters) of empty sites. This is true for p — 1 also since P(0) = in 
this case. Using P(0, 0) = in eq. © and then approximating the rest of the joint distributions 
by the product of single site distributions, we find that while this theory correctly predicts 
the critical density to be one for all p > 0, it gives 5* oc p — p c . Thus, the simple mean field 
analyses above fail to capture the interesting behavior of the activity and underscores the 
importance of correlations in the system. 

Relaxation dynamics. We now discuss the dynamics of the relaxation to the steady 
state for p < p c . At the critical point, we measured the distribution P(t,L,p) of the time 
t required for the system to reach the steady state starting from an initial condition with a 
single empty site. As shown in the inset of fig. the normalised distribution P(t,L,p) is of 
the scaling form 

P(t,L,p)^L- z f(t/L z ) , (10) 

where the scaling function f(x) decays as x~ T for x <C 1 and exponentially for x ^> 1. Since 
the initial configuration can reach the critical state in one time step with probability p/2 
if the empty site is adjacent to the site with two particles, and such an initial state occurs 
with probability 1/L, it follows that P(l, L.p) ~ p/L. Using this in ea. (jl()fl . we find that the 
exponents z and r are related by the scaling relation z(l — r) = 1 with r < 1. Our simulations 
indicate that z ~ 1 + p which matches with the results quoted earlier for the limits p = and 
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Fig. 4 - Semilog plot to show the stretched exponential decay of the activity S(t,p,p) at p — 1/2 for 
p = 1/4(A), l/2(x) and 3/4(+) with L = 2 20 . Inset: Plot of -t~ a In S(t, p,p) vs. i to determine a 
as explained in the text. 



1, and gives r fap/(l + p). As shown in fig. EI the power law decay of P(t,L,p) is consistent 
with this value of r. The above conjectured expression for z was also used to obtain a data 
collapse for the activity to the following scaling form 

S(t,p C:P )^t- 9 g(t/L z ) , (11) 

with p-dependent exponent 9 (see inset of fig. 

We also studied the decay of the activity in the inactive phase with initially Poisson dis- 
tributed particles. Similar to other stochastic sandpiles [18], the activity decays as a stretched 
exponential and as shown in fig.0] we obtain a best fit to the form S(t, p,p) = So exp(— (t/to) a ) 
with a w 0.48,0.42 and 0.36 for p = 1/4, 1/2 and 3/4 respectively. In the inset of fig. El we 
plot the equation 

-ln5 _ - In 5b 1 
t a t a t% ' 

where a is determined by requiring that at large times, the left hand side approach a constant. 
The top curve for p = 1/4 uses a = 1/3 which we recall is the value of the exponent a at 
p = 1. Since this is an increasing function at large times, the true value of a is larger than 
1/3 and in fact, a constant can be obtained with a « 0.394. We take this as evidence to 
rule out a = 1/3 for p < 1. Further, since the bottom curve for p = 3/4 with a « 0.394 
is monotonically decreasing, we conclude that the exponent a is a decreasing function of p 
which agrees with the trend obtained by the best fit above. Thus, for p < p c , the approach 
to the steady state is fast for small p and slows down as p increases. 

We close this discussion by noting that while at p = 1, the relaxation dynamics of the 
activity could be understood via a mapping to the well known two-species annihilation problem 
A + B — > [15], the corresponding reaction-diffusion system for p < 1 turns out to be more 
complicated involving additional reactions such as the ones in which both A and B can be 
created. 
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Conclusions. - We have checked that the critical behavior of our model is robust against 
perturbations. We studied this model with parallel update both with and without multiple 
topplings and found the same critical behavior as for the random sequential case discussed here. 
Further, a preliminary numerical study of the SOC state of our model supports the scaling 
form eq. 1)10(1 for the relaxation time distribution P(t, L, p) with the conjectured z. To conclude, 
we have studied several one-dimensional sandpiles in which the number of particles transferred 
on relaxing an active site is determined stochastically and identified a novel universality class 
of active-absorbing transitions characterised by nonuniversal critical exponents. An analytical 
understanding of this nonuniversal behavior is clearly very desirable. A promising direction 
seems to be a study of the SOC state using the operator algebra introduced in [13]. While the 
largest eigenvalue of these operators correspond to the steady state, the next eigenvalue gives 
the dependence of the relaxation time on the system size L. For our model, we expect a In L 
correction to the usual power law dependence of the next eigenvalue. To see this logarithmic 
function, a study of sufficiently large L is required which would be a subject of future work. 
It may also be interesting to see if this non-universality holds in higher dimensions. 
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